/*
 * MersenneTwister.h
 *
 *  Created on: 31/08/2012
 *      Author: Emeson J. S. Pereira
 */

#ifndef MERSENNETWISTER_H_
#define MERSENNETWISTER_H_



#endif /* MERSENNETWISTER_H_ */

#include <stdio.h>
#include <time.h>
#include <math.h>

//===========Period Parameters=============
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df
#define UPPER_MASK 0x80000000
#define LOWER_MASK 0x7fffffff

//=========Tempering Parameters============
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)

//=========================================
#define NV_MAGICCONST 4 * exp(-0.5)/sqrt(2.0)
//=========================================
double random_uniform( register double a, register double b);
int random_list(int list[], register int size_list);
double normal_variate(double mu, double sigma);

//=========================================
static unsigned long mt[N];
static int mti=N+1;

void
sgenrand(seed)
	unsigned long seed;
{
	mt[0]= seed & 0xffffffff;
	for(mti=1; mti<N; mti++)
		mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
}

double
genrand()
{
	unsigned long y;
	static unsigned long mag01[2]={0x0, MATRIX_A};

	if (mti >= N){
		int kk;

		if (mti == N+1)
			sgenrand(time(NULL));

		for (kk=0;kk<N-M;kk++){
			y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
			mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
		}
		for (;kk<N-1;kk++){
			y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
			mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
		}
		y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
		mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];

		mti = 0;
	}

	y = mt[mti++];
	y ^= TEMPERING_SHIFT_U(y);
	y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
	y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
	y ^= TEMPERING_SHIFT_L(y);

	return ( (double)y / (unsigned long)0xffffffff );
}

double //Get a random number in the range [a, b) or [a, b] depending on rounding.
random_uniform( register double a, register double b)
{
	return a + (b-a)*genrand();
}

int //Choose a random element from a non-empty sequence.
random_list(int list[], register int size_list)
{
	return list[(int)( genrand() * size_list )];
}

double
normal_variate(double mu, double sigma)
{
	/* mu = mean, sigma = standard deviation

       Uses Kinderman and Monahan method. Reference: Kinderman,
       A.J. and Monahan, J.F., "Computer generation of random
       variables using the ratio of uniform deviates", ACM Trans
       Math Software, 3, (1977), pp257-260.
    */
	double u1, u2, z, zz;

	while(1){
		u1 = genrand();
		u2 = 1.0 - genrand();
		z = NV_MAGICCONST*(u1 -0.5)/u2;
		zz = z*z/4.0;
		if(zz <= -log(u2)) break;
	}
	return mu + z*sigma;
}
